Non-equilibrium ribbon model of scroll waves with twist 



Bias Echebarria, 1 Vincent Hakim, 2 and Herve Henry 3 

Departament de Fisica Aplicada, Universitat Politecnica de Catalunya, Doctor Maranon 44> E-08028 Barcelona, Spain. 

Laboratoire de Physique Statistique CNRS-UMR8550, 
Ecole Normale Superieure, 24 rue Lhomond 75231 Paris, France. 
3 Laboratoire de Physique de la Matiere Condensee CNRS-UMR7643, Ecole Polytechnique Palaiseau, France. 

(Dated: February 8, 2008) 

We formulate a reduced model to analyze the motion of the core of a twisted scroll wave. The 
model is first shown to provide a simple description of the onset and nonlinear evolution of the 
helical state appearing in the sproing bifurcation of scroll waves. It then serves to examine the 
experimentally studied case of a medium with spatially-varying excitability. The model shows the 
role of sproing in this more complex setting and highlights the differences between the convective 
and absolute sproing instabilities. 

PACS numbers: 89.75.Kd, 82.40.Ck, 82.40.Bj, 87.19. Hh 



Scroll waves, spirals tridimensional (3D) counterparts, 
are essential structuring elements of the dynamics of 
thick excitable media and are thought to play an impor- 
tant role in ventricular fibrillation^, 0. This has mo- 
tivated detailed examinations of their instabilities, both 
with chemical reactions in gels 0, 0] and theoretically 
HEHESOii- Winfree et al. discovered that twist 
can destabilize a scroll straight core and lead it to adopt 
a helical shape .5|. This "sproing" bifurcation resem- 
bles the twist- induced instabilities of elastic rods [ill 
or of DNA ^jjj , but has remained somewhat of a puz- 
zle since, using long- wave expansions the dynamics 
of a scroll core filament was found to be independent of 
twist 0. Here, insights gained from a numerical stabil- 
ity analysis || lead us to formulate a simple model of 
the core dynamics of a twisted scroll wave that is an- 
alytically tractable and agrees semi-quantitatively with 
the results of reaction-diffusion (RD) simulations. We 
first show that the model provides an easy understand- 
ing as well as an accurate description of sproing. Sys- 
tematic variations of electrophysiological properties are 
known to exist in the heart and gradients of excitabil- 
ity have been shown to promote scroll wave instabilities 
in chemical media |3j.|jJ. Therefore, we then choose the 
case of a medium with spatially varying excitability to 
test the usefulness of the approach beyond the simplest 
case of a homogeneous medium. The results show that 
the observed instabilities 0, are tightly linked to spro- 
ing and illustrate the subtleties brought by the problem 
non-equilibrium setting. 

The center of rotation of a planar spiral becomes for 
a three-dimensional scroll wave the line of instantaneous 
center of rotations R(er, t) where a parameterizes this 
center filament and t is time. We choose to describe the 
scroll wave core as the ribbon (R(c, t), p(tr, t)), where the 
unit vectors p(<x, t) point orthogonally from the center 
filament to the line of instantaneous scroll wave tips. The 
local rotation of the line of instantaneous tips around the 



center filament defines the scroll wave twist t u 



(p x d s p) • T, 



(1) 



where T = dH/ds is the local tangent to the center fil- 
ament and s its arclength (with d s simply a notation 
for l/\dH/da\ d a ). Starting from a straight, untwisted 
scroll wave, a gradient expansion |a, tl la LLj based on 
the translational and rotational invariance neutral modes 
shows that, at lowest order, the scroll core motion_is_sim- 
ply driven by its center filament curvature k 



Ri • N — a\ k, R t • B = a.2 k, 



(2) 



where N is the filament normal, kN = dT/ds, and 
B = T x N its binormal. The case when a\ < is 
analogous to the filament having a negative line tension, 
and the allied instability has been extensively studied 
S El El However, Eq. © leaves twist-induced in- 
stabilities unexplained, since the motion of the mean fil- 
ament is not influenced by the ribbon twist J\. Twist 
appears at higher orders in the gradient expansion of 
ref. 0, Q, @ but, besides being somewhat cumbersome, 
this rigorous approach suffers from the fundamental diffi- 
culty that sproing sets in at a finite wavenumber || . Con- 
sequently, it cannot be precisely described by a gradient 
expansion cut to any finite order [l5| . Therefore, we find 
it instructive to formulate here a simple phenomenolog- 
ical model that captures the essence of the phenomenon 
and that contains only terms essential for the instability 
description. The filament velocity in its normal plane is 
written as a generalization of Eq. @ 

[R-i]_L = fllR-ss + x R ss + diT w R s X R sss (3) 

— d*2T w {Tl sss ]± — biT w [R ssss ]± — 62R.S x Rssss, 

where the brackets denotes the component of the en- 
closed vector orthogonal to the filament tangent (e.g. 
[Rt]± = R t (Rt • T)T). It is worth remarking that 
the helical instability of an elastic ribbon with a gra- 
dient dynamics based on extension and curvature and 



2 



twist energies [l2( essentially depends on the ai,6i and 
d± terms. The 02, 62 and di terms describe motion in the 
orthogonal direction. They can appear in the present 
non-potential problem due to the handedness of the spi- 
ral rotation and their sign depends on the spiral sense of 
rotation. Eq. © needs to be completed by the evolution 
of the ribbon twist. The twist kinematics can be adapted 
from previous investigations of elastic ribbons. Following 
ref. [la], we note that as one slides along the central fil- 
ament at a fixed time t, the ribbon vector p rotates and 
remains orthogonal to the filament tangent T, 



dp 
~da~ 



ds 
da 



(T x p) - 



dT 
~da~ 



p T. 



(4) 



This is also true as time evolves when one stands at a 
fixed abscissa a and similarly, 



dp 

~di 



-£=u(Txp)- — p T 



dT 

di 



(5) 



where lu is the local instantaneous scroll wave rotation 
frequency. Comparing crossed-derivatives of Eqs. |0J 0) 
gives as single compatibility condition the equality of the 
projections of dt.aP and d G: tP on T x p, 



d_ 

dt 



ds 

da 



duj 
~da~ 



dT dT 

~da X ~dt 



(6) 



The kinematic Eq. © is a local description for an exten- 
sible ribbon of the well-known conversion of twist into 
writhe 01 associated to linking number conservation at 
a global level. The specific dynamics of the present prob- 
lem is encoded in the twist-dependent rotation frequency 
lu . A good approximation for moderate twist is obtained 
by keeping the first twist corrections to the untwisted 
scroll frequency loq, 



w 



Dd s 



(T • 5 t R)r. u 



(7) 



where the coefficients D and c can be explicitly calculated 
by linearization around the straight scroll wave and pro- 
jection over the adjoint eigenmodes [3. The last term in 
Eq. J7J) is due to the apparent rotation of p coming from 
changing position along the filament. Eq. JHJ) with J7|) is 
equivalent to Keener's phase equation [6| and completes 
our formulation of the ribbon model. In the following, 
this simplified model is compared to simulations of RD 
equations in the form 

d t u — V 2 u + u(l-u)[u-(v + j3)/a]/e, d t v = u-v, (8) 

with a = 0.8, e = 0.025, and different values of j3. 
Eqs. (jHJ are simulated with an explicit second order 
scheme, with dx = 0.15, and dt = 5.625 • 10 -3 . 

- Sproing. Taking a vertical filament along the z-axis 
and assuming small transverse X, Y deformations, Eqs. 
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FIG. 1: a) Dispersion relations obtained from a direct linear 
stability analysis of a twisted scroll wave of Eq. © H, with 
a = 0.8, e = 0.025, and = 0.01, and b) from Eq. @, 
with a — 0.2 + 0.2i, d — 3.5 + i, and b — 2 + i, chosen to 
give a semi-quantitative overall agreement between the two 
sets of curves. Above a crit ical twist, Re(S7) is positive for k 
around a non-zero k c and the straight scroll becomes unstable 
to a finite wavenumber Hopf instability. Since the instability 
occurs at finite wavelength, a better fit at k = (using the 
exact value of a ||) typically deteriorates the overall fit. 



(I.'')l(il7ll become to quadratic order, 

d t W = ad 2 z W + idr w d 3 z W - bd^W, (9) 
d t r w = d s (Dd s r w ) + d s (cr 2 ) + d s uj 

+Rc{[d z {T w d z W) - i(a%W)d,]dtW}, (10) 

where a complex notation has been used for the deforma- 
tion field W(z,t) — X(z,t) + iY(z,t) and the constants 
a,b,d (e.g. a = a% + ia%), and d s can be approximated 
by [1 — \d z W\ 2 /2]d z . For a uniformly twisted filament 
in a homogeneous medium (loq = est.), the linear modes 
W(z, t) — e lfez+0t , correspond to helices of pitch k. Their 
dispersion relation is obtained from Eq. © as, 



VL = -ak 2 + dr w k 3 -6fc 4 . 



(11) 



With appropriate constants a,b,d, it is similar to the dis- 
persion relation obtained from RD Eq. (JSJ (Figs. Q^,b). 
A secondary local maximum appears away from k = 
when t w = A/Zy/2a\b\/d\. In a large box, instability 
sets in at the critical twist r£, = 2y / aibi/d 2 , with the 
pitch of the allied helix equal to k c — a\jb\. The mode 
k becomes unstable above tJ 
a homogeneous twist v, 
R(t) of a helix of pitch k grows as 



„. k when Re[n(k)} > 0. For 

w,k' 



slightly above rf„ u , the radius 



Rt = lk(r w - T^ k )R, 



(12) 



where 7^ = Re[dfl(k) / dr w ] = d\k 3 . Saturation of the 
instability comes from the coupling between twist and 
bending described by Eq. (|10fl . For an helical mode of 
pitch k and time dependent but homogeneous twist and 
radius, the partial s derivative terms of Eq. i|10fl vanish. 
The last and only remaining term is equal to —k (r w + 
k)RR t so integration of Eq. (|10(l shows that the twist t w 
decreases with the helix radius, 



(r»+k)k 2 R 2 /2, 



(13) 



FIG. 2: Distributions of twist along a straight filament in 
simulations of RD Eq. J8J (solid line), and given by Eq. dTUJ 
(dashed line), a) Jump of excitability obtained by taking 
@{z) = P b + (fa - Pb)Q{z - L/2) in Eq. © or in the model 
by taking the allied spiral frequency jump u>o = uit + {^t — 
uob)Q{z — L/2). b) Linear gradient of excitability: /3(z) = 
fib + {fit - f3b)z/L in Eq. © or uj = u b + (w< - u>b)z/L for 
the model. The parameters are the same as in Fig. Q with 
b =O.Ol, fit = 0.03, Lo b = 1.80, u t = 1.696, and coefficients 
D = 0.578, c = 0.720, at the bottom, and D = 0.614, and c = 
0.856, with equivalent expressions. This gives a theoretical 
prediction of r™ ax ~ (Aw /c) 1/2 = 0.35 in a). Here and in 
the following we plot twist i n absolute value. 

where r° is the initial twist of the straight scroll, and we 
have assumed (Rk) 2 <C 1 0|. Comparison of Eq. I|12(l 
and ljT3|l describes sproing as a supercritical bifurcation 

Rt = J k (r° - T° >k )R -^(r°+ k)k 2 R 3 . (14) 

The deformation of the center filament decreases the ini- 
tial twist until the critical value t w = k , is reached 
at which point the driving force for the instability disap- 
pears. The final helix radius is R — [2(r° — fc )/(r° + 
/c)^ 2 ] 1 / 2 (with k — k c in a large box). 

These analytic results compare well to results of RD 
simulations with periodic boundary conditions (BC) in 
the z-direction to enforce linking number conservation 
[20|. As previously reported sproing is found to be 
a supercritical bifurcation and the twist of a bifurcated 
helix is very close to the critical one, in good agreement 
with the above findings. In large systems, as for oscilla- 
tory media .2.1 J , the helices resulting from sproing may be 
unstable to secondary Hopf instabilities || which appear 
sensitive to higher order nonlinearities not included in 
Eq. J2J). These can be described by amplitude equations 
for the coupled helix amplitude and excess local linking 
number. The equations can be derived from the reduced 
model or directly from the RD Eq. (|5J) and take a form 
similar to other cases with a conservation law [j^. In 
simulations of Eq. © these secondary instabilities typi- 
cally result in other helices with smaller wavenumber, or 
in modulated structures. 

- Inhomogeneous twist. Most experimental situations 
correspond to imposing free non-flux BC on Eq. (JHJ) 
rather than periodic ones. These do not conserve to- 
tal linking number and an initially twisted scroll wave 
untwists |3( in an homogeneous medium. Spatial varia- 
tions of excitability do however promote twist formation. 



FIG. 3: a) Center filament maximal radius vs fi t . The straight 
scroll becomes unstable at fit ~ 0.03939, resulting in a small 
hysteresis, b) A 3D view of the solution for fit = 0.0399. 



Fig. [2] shows RD simulations for a straight vertical fil- 
ament, with non-flux BC, and either a jump (Fig. 
or a linear gradient (Fig. 03) of the value of fi in the 
z-direction. For moderate variations of fi, the initial un- 
twisted core, remains straight but evolves toward a final 
twisted and steadily rotating configuration. The different 
natural spiral rotation frequencies ujq(z) create phase dif- 
ferences between different heights z which together with 
the rotation frequency increase with twist [Eq. Q] lead 
to a steady state. The resulting distribution of twist can 
be computed, either analytically or numerically, from the 
model Eq. I|10|) using the appropriate source term d z ujQ. 
The calculated distribution of twist agrees remarkably 
well with the RD simulations as shown in Fig. For 
definitcness, we focus on the case of a medium with an 
excitability jump. In the RD Eqs. (JSJ), we fix fib = 0.01 
in the medium bottom half part, and take f3 — fit > (3b in 
the less excitable top half. When the jump in excitability 
is larger than a critical value, the straight scroll becomes 
unstable, ve ry s imilarly to what is observed in experi- 
ments U 0, [13 ■ The instability is slightly subcritical 
and the resulting structures modulated helices (Fig. |3J|. 
To clearly relate this instability to sproing, we consider 
now the limit of large systems. Then, for a moderate 
jump of excitability, the scroll core is straight and its 
frequency is basically set by the domain most excitable 
half where the scroll twist is negligible. The scroll twist 
T max j n ^ e domain less excitable half is almost constant 
and simply determined by the frequency jump Aw be- 
tween the two domain parts, r™ ax ~ (Auiq/c) 1 ^ 2 . When 
T jnax reacnes t ne sproing threshold, for a large enough 
jump, one could expect sproing to set in with the center 
filament taking the shape of a helix of constant radius 
in the low excitability region and radius decreasing to 
zero in the higher excitability part. However, the insta- 
bility onset differs from the sproing threshold in a ho- 
mogeneous system, even when L — > oo (Fig. ^i). Fur- 
thermore, the bifurcated filament radius decreases expo- 
nentially also in the region of constant twist (Fig. 
In order to clarify the phenomenon, we have analyzed 
the ribbon model in this geometry. We have solved the 
eigenvalue problem given by Eq. @, with the distri- 
bution of twist calculated with Eq. i|10[l (with constant 
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FIG. 4: a) Instability onset vs domain size L, for domains with 
a jump in excitability. For fixed L, (3t is increased until the 
system becomes unstable. The corresponding critical twist 
(filled circles) is plotted together with the fit (dashed line): 
T max = 0. 344 + 31. 3/L 2 . For comparison, the onset of sproing 
for periodic BC (r^ ~ 0.29) is also shown (dotte d line) for 
13 = 0.03, that corresponds to the critical value of f3t in a lar ge 
system, b) Solutions for f% = 0.01 and (3 t = 0.031. c) Model 
critical twist obtained by solving the eigenvalue problem vs L 
for the distribution of twist shown in Fig. [3^ (diamonds), a 
constant twist in half of the system (squares), and vs L/2 for 
a constant twist in the whole system (circles). For comparison 
we also show the onset of convective (dotted line), absolute 
instability (long dashed line), and the onset obtained with the 
double root criterium (short dashed line), d) Critical modes 
obtained from the eigenvalue problem, using the distribution 
of twist, for L = 150. 



values of D = 0.578 and c = 0.72). Similarly to RD 
simulations, an instability develops in the region of con- 
stant twist when r™ ax is large enough but its threshold 
differs from the sproing threshold for periodic domains 
(Fig. 13k). The instability is nonetheless related to spro- 
ing. The reason is that periodic BC allow the growth 
of convective instabilities, that decay with non-flux BC. 
For a complex growth rate Q, four complex wavenum- 
bers fcj(Q) satisfy the dispersion relation Eq. illlfl . The 
relevant sproing absolute spectrum, for a given constant 
twist in a large domain, lies on the curve of complex 
such that Im[k2(£i)] = Jm[fca(f2)l, with the hi ordered 
by increasing imaginary part |24ll25|. For low twist, this 
curve lies entirely in the Re(Q) < half plane. The 
absolute instability threshold twist, for which the curve 
crosses the fl imaginary axis [2(| , coincides with the large 
L limit of t c (L) as shown in Fig. QJ:. The most unstable 
modes at threshold are two counterpropagating waves, 
with the same spatial growth rate, and nonzero group 
velocity. The similarity between the critical modes for 
the ribbon model and the RD equation (Fig. fur- 
ther shows that sproing is also a likely explanation of the 
latter case and of the experimental observations 0, 0] . 
In conclusion, the proposed ribbon model provides a 



semi-quantitative description of the motion of twisted 
scroll waves and a clear understanding of several features 
that are difficult to directly extract from RD equations. 
This will hopefully help to further analyze scroll wave dy- 
namics in complex media and to better assess the effects 
of gradients of electrophysiological properties and other 
complicating features in the cardiac muscle. 
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FIG. 1: Simulations of reaction-diffusion Eq. (8) of the main 
text. The core of the initial twisted scroll is straight and 
vertical, and periodic boundary conditions are applied in the 
z-direction to enforce linking number conservation, c) Resta- 
bilized scroll mean filament, for the parameters correspond- 
ing to Fig. 1 a) of the main text (a — 0.8, e = 0.025, and 
/3 = 0.01), and t w = 0.45. d) Nondimensional radius of the 
helix as a function of the reduced distance to the onset of 
sproing, for several values of /3 (full symbols) in a small box 
of height L with a single turn of helix (r^ = k — 2-k/L). 
The predicted straight line of slope one (dotted) is shown for 
comparison. Note that this linear relation holds beyond small 
departures from threshold since linking number conservation 
and restabilization of the bifurcated helix at the critical twist 
give 

T^h = t°/[1 + (Rt°) 2 ] without any expansion, as follows 
from footnote [19] of the main text with k — t%. In the limit 
(Rt°) 2 < 1 this expression reduces to the one obtained from 
Eq. (13). 



